# Reconstructed visualization code for Figure 6B-C.
# Uses exported result tables, NOT raw read counts; does not rerun edgeR.
# Run: Rscript Supplemental_File_1_Figure6BC.R [bundle_directory]
# Base R only. Original historical plotting settings were not supplied.
args <- commandArgs(trailingOnly = TRUE)
script_arg <- grep('^--file=', commandArgs(), value = TRUE)
root <- if (length(args)) args[1] else if (length(script_arg)) {
  dirname(normalizePath(sub('^--file=', '', script_arg[1])))
} else getwd()
out <- file.path(root, 'plots')
dir.create(out, recursive = TRUE, showWarnings = FALSE)
selected <- read.csv(file.path(root, 'data', 'Figure6B_selected_piRNAs.csv'),
                     check.names = FALSE, stringsAsFactors = FALSE)
all <- read.csv(file.path(root, 'data', 'Figure6C_all_piRNAs.csv'),
                check.names = FALSE, stringsAsFactors = FALSE)
samples <- c(paste0('Myopia', 1:5), paste0('Control', 1:5))
required <- c('ID', 'logFC', 'pvalue', 'qvalue', samples)
for (d in list(selected, all)) {
  stopifnot(all(required %in% names(d)), !anyDuplicated(d$ID))
  stopifnot(all(is.finite(as.matrix(d[, c('logFC','pvalue','qvalue',samples)]))))
  stopifnot(all(d$qvalue >= 0 & d$qvalue <= 1), all(d$pvalue >= 0 & d$pvalue <= 1))
  stopifnot(all(as.matrix(d[,samples]) >= 0))
}
q_cutoff <- 0.05
fc_cutoff <- 1
sig <- all$qvalue < q_cutoff & abs(all$logFC) > fc_cutoff
stopifnot(setequal(selected$ID, all$ID[sig]))
idx <- match(selected$ID, all$ID)
stopifnot(isTRUE(all.equal(selected$logFC, all$logFC[idx])),
          isTRUE(all.equal(selected$qvalue, all$qvalue[idx])))
# The values supplied in sample columns are normalized expression values.
# New, explicit visualization choice: log2(expression+1), then row z-score.
# Set transform <- 'none' for row scaling of untransformed expression.
transform <- 'log2p1'
mat <- as.matrix(selected[,samples]); rownames(mat) <- selected$ID
if (transform == 'log2p1') mat <- log2(mat + 1)
sd_rows <- apply(mat, 1, sd)
if (any(sd_rows == 0)) stop('Constant-expression rows cannot be z-scored.')
z <- t(scale(t(mat)))
hc <- hclust(dist(z, method='euclidean'), method='complete')
ordered <- z[hc$order,,drop=FALSE]
write.csv(ordered, file.path(out,'Figure6B_row_zscores.csv'))
pal <- colorRampPalette(c('#234BA2','white','#EE2424'))(101)
lim <- max(abs(z))
pdf(file.path(out,'Figure6B_heatmap.pdf'), width=8, height=7, useDingbats=FALSE)
layout(matrix(c(0,1,5,2,3,4),nrow=2,byrow=TRUE),widths=c(1.1,5,1.3),heights=c(.7,6))
par(mar=c(0,0,2,0))
plot.new();plot.window(xlim=c(.5,10.5),ylim=c(0,1),xaxs='i',yaxs='i')
rect((1:10)-.5,0,(1:10)+.5,1,col=c(rep('#F18C85',5),rep('#42B9C1',5)),border=NA)
title(main='Heatmap',line=0)
par(mar=c(5,0,0,0))
plot(as.dendrogram(hc),horiz=TRUE,axes=FALSE,leaflab='none',
     xlim=c(max(hc$height),0),ylim=c(.5,nrow(z)+.5),yaxs='i')
par(mar=c(5,0,0,0))
image(1:10,1:nrow(z),t(ordered),col=pal,zlim=c(-lim,lim),axes=FALSE,
      xlab='',ylab='',xlim=c(.5,10.5),ylim=c(.5,nrow(z)+.5),xaxs='i',yaxs='i')
axis(1,at=1:10,labels=samples,las=2,tick=FALSE,cex.axis=.8)
par(mar=c(5,1,0,3))
plot.new();plot.window(xlim=c(0,1),ylim=c(-lim,lim))
br <- seq(-lim,lim,length.out=102)
rect(0,head(br,-1),.25,tail(br,-1),col=pal,border=NA)
axis(4,at=pretty(c(-lim,lim)),las=1,cex.axis=.7)
mtext('Row z-score',side=4,line=2,cex=.7)
par(mar=c(0,0,0,0));plot.new()
legend('center',legend=c('Myopia','Control'),fill=c('#F18C85','#42B9C1'),bty='n',cex=.7)
dev.off()
# Preserve supplied logFC sign; never reverse it to force a biological claim.
# Positive/negative refer to the stored model contrast. Confirm the historical
# design/contrast before calling these Myopia-up or Control-up in the manuscript.
status <- ifelse(sig & all$logFC > 0,'Positive logFC',
                 ifelse(sig & all$logFC < 0,'Negative logFC','Not significant'))
levels_status <- c('Positive logFC','Negative logFC','Not significant')
counts <- table(factor(status,levels=levels_status))
cols <- c('Positive logFC'='#F42535','Negative logFC'='#2856A3','Not significant'='#BDBDBD')
y <- -log10(pmax(all$qvalue,.Machine$double.xmin))
pdf(file.path(out,'Figure6C_volcano.pdf'),width=7,height=5.5,useDingbats=FALSE)
par(mar=c(4.5,4.5,2.5,8),xpd=FALSE)
plot(all$logFC,y,pch=16,cex=.5,col=cols[status],
     xlab=expression(log[2](fold~change)),ylab=expression(-log[10](q-value)),
     main='Volcano plot of piRNAs')
abline(v=c(-fc_cutoff,fc_cutoff),h=-log10(q_cutoff),lty=4)
legend(par('usr')[2]+.1,par('usr')[4],xjust=0,xpd=NA,bty='n',cex=.65,pch=16,
       col=cols,legend=paste0(levels_status,': ',as.integer(counts)))
dev.off()
all$plot_class <- status
write.csv(all,file.path(out,'Figure6C_plot_data.csv'),row.names=FALSE)
write.csv(data.frame(class=levels_status,n=as.integer(counts)),
          file.path(out,'Figure6C_counts.csv'),row.names=FALSE)
writeLines(c(paste('Heatmap transformation:',transform),
             'Row z-score; Euclidean distance; complete linkage; fixed sample order.',
             'Thresholds: qvalue < 0.05 and abs(logFC) > 1.',
             'Model contrast direction: retained as supplied; requires original design records.',
             paste('Input MD5:',names(tools::md5sum(list.files(file.path(root,'data'),full.names=TRUE))),
                   tools::md5sum(list.files(file.path(root,'data'),full.names=TRUE)))),
           file.path(out,'parameters.txt'))
capture.output(sessionInfo(),file=file.path(out,'sessionInfo.txt'))
print(counts)